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ABSTRACT 



Aims. We have recently developed a microscopic Monte Carlo approach to study surface chemistry on interstellar grains 
and the morphology of ice mantles. The method is designed to eliminate the problems inherent in the rate-equation 
formalism to surface chemistry. Here we report the first use of this method in a chemical model of cold interstellar cloud 
cores that includes both gas-phase and surface chemistry. The surface chemical network consists of a small number of 
diffusive reactions that can produce molecular oxygen, water, carbon dioxide, formaldehyde, methanol and assorted 
radicals. 

Methods. The simulation is started by running a gas-phase model including accretion onto grains but no surface 
chemistry or evaporation. The starting surface consists of either fiat or rough olivine. We introduce the surface chemistry 
of the three species H, O and CO in an iterative manner using our stochastic technique. Under the conditions of the 
simulation, only atomic hydrogen can evaporate to a significant extent. Although it has little effect on other gas-phase 
species, the evaporation of atomic hydrogen changes its gas-phase abundance, which in turn changes the flux of atomic 
hydrogen onto grains. The effect on the surface chemistry is treated until convergence occurs. We neglect all non-thermal 
desorptive processes. 

Results. We determine the mantle abundances of assorted molecules as a function of time through 2 x 10"" yr. Our 
method also allows determination of the abundance of each molecule in specific monolayers. The mantle results can be 
compared with observations of water, carbon dioxide, carbon monoxide, and methanol ices in the sources W33A and 
Elias 16. Other than a slight underproduction of mantle CO, our results are in very good agreement with observations. 

Key words. ISM: abundances - ISM: molecules-molecular processes 



1. Introduction 

Not only are interstellar grain surfaces believed to be 
the sites where molecular hydrogen is formed in diffuse 
clouds, but in cold cores of dense clouds interstellar grain 
chemistry is an important mechanism for the formation of 
quite a few additional species, most of which are produced 
through gra dual hydrogenat ion of gas-phase molecules 
(iTielens fc H agcn 1982; Hasc gawa et al.lll992l:lHerbstlll995t 
lEhrenfreund fc Schutte.200Q) . Methanol is formed in such 
a manner, and its surface formation has received much 
attention over the years, b o th fr o m theoreticians and 
model ers (iTielens fc Whittetl [l996l: IStantcheva fc Herbstl 



12004 lAwad et al.l 120051 [Garrod et al.l I2006D and ex rori 



mentalists (iHiraoka et al.[ 2002t Watanabe fc Kouchlll2002t 
IWatanabe et al.ll2003f ). 

In order to simulate all of the chemistry that occurs in 
a cold region, it is necessary to include both the chemical 
kinetics in the gas and on the surfaces of dust grains. The 
basic method for simulating chemistry is to use coupled 
rate equations, one for each species in the simulation. Rate 
equations have been proven to be an accurate and efficient 
appro ach to simulat e gas phase chemical kinetics (i HerbstI 
Il995t [Roberts et all |2004 TWakelam et al.l [2006). In these 
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equations, the average concentration of a species is pro- 
gagated forward in time through a non-linear differential 
equation that contains both source and sink terms corre- 
sponding to chemical reactions. Their success comes from 
the fact that average concentrations are normally well de- 
fined in the gas phase as functions of time. In other words, 
the so-called mean-field approximation is valid. However, a 
straightforward and efficient generalization to simulate dif- 
fusive su rface chemistry on interst e llar grains in a similar 
manner ("Pickles & Williams 197'/: 'Hascga wa et al.l Il992f ) 
can be in error (Tielens & Hagcn 1982). In particular, when 
the number of surface species on an average grain is small, 
the rate equation method breaks down since its mean-field 
approach does not accurately describe the surface kinetics, 
which must be described by a method in which both the dis- 
crete number of species and fiuctuations in this number are 
considered. Several alternative methods to rate equations 
have been proposed to overcome this difficulty. The simplest 
of these is the modified-rate-equation approach, which uses 
the rate equation method but instead of using the actual 
diffusion rate on the surface, scales this rate downward so 
that the rate of reaction does n ot on average exce ed the 
rate of accretion of a reactant ([Caselli et al.l[T998[ l. This 
semi-empirical approach is efficient c omputationally, but it 
is not reliable under all con ditions (jStantcheva fc Herbstl 
12004 ICxarrod fc Herbstl l200fih . 
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There are two more rigorous approaches, both of which 
are stochastic and macroscopic in nature. In other words, 
they allow a determination of the probability that a certain 
number of molecules of a given species exists on a grain at 
any time without distinguishing where on the gra in the ad- 
sorba te might lie. In the Monte Carlo approach (jCharnlevI 
Il998( ). random numbers are used to determine what pro- 
cesses occur in a tim e interval based on their relative rates 
(jCaselli et al.l^002l ). It is difficult to develop a combined 
gas-grain model, in which the Monte Carlo method is used 
for the surface chemistry and the rate-equation approach 
for the g as-phase chemist ry. The master equation appr oach 
(|Biham et al..,200L : , Green at al.l[200ll : ICharnlevll200lf ). on 
the other hand, is based on integrating a set of ordinary 
differential equations for the probabilities of surface pop- 
ulations, and can be coupled with rate equations for gas 
phase chemistry. Indeed, a coupling of the master equation 
approach to surface chemistry and a Monte Carlo approach 
to ga s-phase chemistry has also been suggested (jCharnlevI 
Il998l) . Although the master equation approach is believed 
to be the better candidate to solve the difficulties with the 
surface rate-equation approach, it has the drawback that 
there are many simultaneous equations to be solved, since 
the probabilities of surface populations for different species 
can be correlated. Nevertheless, with suitable approxima- 
tions, two groups have successfully run gas-grain networks 
of cold interstellar cloud cores with the master equation 
approach to model surface che mistry ([Stantcheva fc HerbstI 
|2004 [Lipshtat fc Bihamll2004l ). In these simulations, only a 
small subse t of the surface reactions used in full-size mod- 
els (.Hasegawa et al.lll992t) were included, with methanol 
the most complex species synthesized. 

Recently, we adopted a microscopic stochastic method 
for surface chemistry, known as contin uous-time random 
walk (CTRW) Monte Carlo approach (|Montroll fc Weisd 
Il965() . which tracks the specific trajectory and rate of dif- 
fusion of each adsorbate and the binding (desorption or 
evaporati on) energy and barrier to d iffusion a t each lo- 
cal site (jChang et all l2005t ICuppen fc Herbst, |2005( ). In 
this method, grain surfaces are modeled as square lattices 
while hopping, evaporation and deposition are regarded 
as Poisson processes. We developed the CTRW approach 
for several reasons. First, for a surface with a continu- 
ous distribution of diffusion barriers and evaporation en- 
ergies, neither the rate equation approach nor the macro- 
scopic stochastic approaches can be easily implemented. 
Even for a surface with discrete values of the energy pa- 
rameters, known as a rough surface, the earlier approaches 
are not readily useable. Moreover, it is facile to consider 
time-dependent barriers and desorption energies with the 
CTRW approach. Secondly, all previous methods neglected 
the "back diffusion" of hopping species on a granular sur- 
face. This neglect, however, turned out to be relatively 
unirnportant for astron omical purposes (|Lohmar fc Krua 
[20061 lChang"eraIll200l . Finally, the use of a microscopic 
Monte Carlo method allows us to follow the morphology of 
a mantle as it develops, so that only reactions in nearby 
monolayers can react with one another. For instance, a hy- 
drogen atom on the surface is not likely to react with an 
oxygen atom deeply buried in the ice mantle. This con- 
straint has been neglected in almost all previous gas-grain 
network calculations. 

Although the CTRW Monte Carlo approach allows us 
in principle to study both the detailed microscopic chem- 



istry occurring in a grain mantle and the actual mor- 
phology of the mantle, use of the method presents some 
formidable computational challenges. Up to the present, we 
have used the method only to study the details of molecular 
hydrogen formation on surfaces releva. i it to diffuse clouds 
(IChang et al.l l2005t ICuppen fc HerbstI 120051 IChang et al.l 
120061). In a work in preparation, we are expanding our use 
of the method to study the formation of mantles of water 
ice and related species whe n, hydrogen and oxygen atoms 
accrete onto a grain surface (jCuppen fc HerbstlfeoOTf ) . This 
extension still does not represent a true gas-grain model 
because the gas-phase chemistry is not followed. Indeed, 
there is no straightforward way to implement varying gas- 
phase abundances at the same time the surface chemistry 
progresses, which makes it difficult to use the CTRW ap- 
proach simultaneously with a gas phase kinetics simulation 
program. 

The difficulty with time-varying gas-phase abundances 
is that these lead to time- varying fluxes of accreting species 
onto the grains. The CTRW method calculates deposition 
as a Poisson process with a time interval At between events 
calculated from a random number X between and 1 de- 
termined via the relation At = —tlnX, where t is the aver- 
age time interval between deposition. If abundances change, 
this equation would be invalid since t would no longer be 
constant and the deposition would not be a strictly defined 
Poisso n process. This problem has been solved bv I Janse"nl 
(|1995D for the simulation of a laboratory procedure known 
as temperature-programmed desorption (TPD), in which 
the temperature changes, resulting in a change of residence 
time for surface species. This technique, however, cannot be 
directly applied to a gas-grain chemical model because, un- 
like the TPD simulation, where the temperature is linearly 
dependent on time, no prior knowledge of the abundances 
of gas-phase species is available. 

In the present paper, we present a first attempt to com- 
bine a treatment of surface chemistry by the CTRW ap- 
proach with a treatment of gas-phase chemistry by the 
standard rate equation approach. The iterative method is 
based on rather weak coupling between the two chemistries, 
a condition that occurs primarily for cold cores, where 
thermal evaporation is negligible for most surface species. 
The present paper uses a variant of the so-ca lled H, O 
and CO system as the surface re action network (jCharnlevI 
I2OOII: IStantcheva fc HerbstI [200I . In this system, the gas- 
phase species atomic hydrogen, atomic oxygen, and car- 
bon monoxide accrete onto grains and react to form species 
such as molecular oxygen, water, carbon dioxide, formalde- 
hyde, methanol, and assorted radicals. The surface net- 
work, which is complex enough to mimic a full surface 
chemical reaction network, has been extensively studied by 
the modified rate equation approach (CascUi et al. 2002), 
the macroscopic Monte Carlo meth od (ICaselli et al.l|2002l) . 



and the master equation approach (jStantcheva et al.ll2002| ) 
with fixed gas-phase abundances of the accreting species. 
The master equation has also been used in a calculation, 
like the present, where ga s -phase chemistry also o ccurs 
(jStantcheva fc HerbstI [200I iLipshtat fc BihamI I2004D . By 
carefully studying the infiuence of the desorption of species 
from the grain surfaces on the overall gas phase abundances, 
we are led to our weak coupling hypothesis that the gas- 
phase kinetics is almost independent of the surface chemical 
kinetics. The weak coupling allows us to handle the chem- 
istry in an iterative manner. 
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The paper is organized as follows. Section 2 reviews the 
H, O and CO system, while Section 3 discusses the initial 
conditions and basic assumptions used in our treatment. 
The details of our method of simulation and surface model 
are discussed in Section 4. We present and discuss our re- 
sults in Section 5, including a comparison with observations 
of molecular ices in cold regions. Our conclusions are stated 
in Section 6. 

2. The H, O, and CO system 

In the stochastic simulations to be discussed later, we fol- 
low only the accretion of H, O and CO from the gas-phase 
and the surface reactions that result from this accretion. 
As will be discussed in detail, the accretion of all species is 
considered in a calculation prior to the Monte Carlo simu- 
lation. 

Table [1] shows the limited surface reaction network 
used in our stochastic simulations. With two exceptions, 
the reactions can occur through the diffusive (Langmuir- 
Hinshelwood) mechanism until the reactants find one an- 
other in the same site. The Eley-Rideal mechanism, in 
which a gas-phase species lands atop an adsorbed one and 
reacts with it, is also considered. Although the formation of 
molecular hydrogen via recombination of hydrogen atoms 
is listed, as we discuss below, this process need not be con- 
sidered if one starts with hydrogen in its molecular form. 
Most of the reactions in Table [1] occur without activation 
energy barriers; for those tha t do possess a c tivati on en- 
ergy, we use values taken from Case lli et ahl (j2002t ). The 
H, C, and C O network in Table [11 differs somewhat from 
that used by IStantcheva fc Herbstj ()2004[ ). These authors 
included two reactions to form carbon dioxide - the associ- 
ation reaction between CO and O and the reaction between 
O and the radical HCO - which are relatively unimpor- 
tant in our calculations mainly because the large activation 
energy cannot be efficiently tunneled under by the heavy 
reactants. We also added a few reactions to make the sys- 
tem more realistic. The most important reactions in this 
network concern the gradual hydrogenation of CO to pro- 
duce methanol, which has been studied in th e laboratory 
(jWatanabe fc Kouchi|[200^ iFuchs et al.ll2007D . 

In our simulations, which pertain to a source at a tem- 
perature of 10 K or 15 K, we assume that only H and O can 
hop, and ignore the movement of all other species. The rate 
of hopping depends exponentially on the energy barrier be- 
tween adsorption sites; this barrier is significantly lower for 
these physisorbed atoms than for the other surface species 
followed. At higher temperatures, the hopping of heav- 
ier sp ecies begins to become important ()Garrod fc Herbsti 
|2006[ ). Likewise, thermal evaporation of almost all species is 
unimportant because it depends exponentially on the des- 
orption energy, and this energy is much too large for species 
heavier than hydrogen at 10-15 K. In fact, we need only fol- 
low the evaporation of hydrogen atoms. Both the hopping 
barrier and the evaporation energy for H and O atoms are 
not constant during our simulation. 

Since we follow the morphology of the growing mantle, 
we find that among monolayers that lie next to one another, 
products of a reaction in one of the monolayers may lie atop 
other reactive species in the next lower monolayer and be 
able to react with them. The OH + H2CO ^ HCO + H2O 
reaction in Table [T] is a good example. The atom O cannot 
react with H2CO, but it reacts with H to form OH. If an 



Table 1. Surface reactions in the H, O and CO sys- 
tem 



Number 


Reaction 


Ea{K) 


1 


H + H ^ Ha 




2 


H + ^ OH 




3 


H + OH -> H2O 




4 


H -1- CO HCO 


2500 


5 


H -1- HCO ^ H2CO 




6 


H + H2CO H3CO 


2500 


7 
8 


H + H3CO ^ CH3OH 
+ ^ O2 




9 


+ OH O2 + H 




10 


CO -f OH CO2 + H 


176 


11 


OH + H2CO -> HCO + H2O 




12 


+ H3CO ^ H2CO + OH 





H2CO molecule lies below the original oxygen atom, the 
newly-formed OH product will directly react with H2CO. 
What is interesting is that OH and H2CO are assumed to 
be stationary under these conditions, but they can react by 
the diffusion of H atoms. 

3. Initial conditions and some basic assumptions 

It is common in pseudo-time-dependent models of dense 
cloud chemistry to assume that the initial form of hydro- 
gen is molecular, having been formed on grain surfaces dur- 
ing the diffuse cloud stage or somewhat later as the dense 
core was formed. In dense cloud cores, H2 is consumed to 
some extent by the gas phase chemistry. Indeed, previous 
model calculations at 10 K show that cosmic ray ionization 
of H2 leads to a residual atomic abundance of atomic hy- 
drogen of ~ 1 cm""^ while that for H2 is maintained at 10'* 
cm^'^. But does this maintenance depend on the continu- 
ing surface formation of H2 under dense cloud conditions? 
Assuming a standard grain size of 0.1 and a standard 
gas-to-dust ratio, a sticking coefficient of unity for H atoms, 
and unit efficiency for the conversion of H to H2 on grain 
surfaces, we obtain that at most 10% of the initial H2 abun- 
dance can be produced within 10^ yr. For a more reasonable 
time scale of 10^ yr, the amount of H2 produced is 100-fold 
less. Despite this minimal production, no depletion of H2 
has ever been reported even at higher temperatures when 
the formation efficiency drops significantly, which indicates 
that gas-phase chemistry during the dense cloud stage does 
not change the abundance of H2 significantly. So, by this 
argument, molecular hydrogen formation on dust grains in 
dense clouds can be neglected, because it does not affect 
the gas. 

In order to confirm the unimportance of H2 forma- 
tion and to understand the strength of the coupling 
between gas-phase chemistry and surface chemistry, we 
performed two calculations at 10 K using the rate equation 
approach. In the first calculation, desig nated Model 
1, we used the osu.2005 gas-grain code (|Garrod et al.l 
2006), which includes over 4000 gas-phas e reac tions (see 
http://www.physics.ohio-state.edu/~eric/research.html) 
and hundreds of surface reactioiisj both of which are 
treated by the rate-equation method. Although all gas- 
phase species were allowed to accrete with unit sticking 
efficiency, we included only those surface reactions in Table 
[1] In addition, atomic hydrogen was allowed to evaporate 
in competition with reaction. In the second calculation. 
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Table 2. Initial Elemental Abundances 



Species 


72i /tIpj 
1 III 1 1 (ij-j^ 


He 


0.09 


o 


3.2 X lO"* 


H 


5.0 X 10"^ 


Ha 


0.5 


c 


1.4 X 10"" 


S 


1.5 X 10"*^ 


Si 


1.95 X 10"® 


Fe 


7.4 X 10"'' 


Na 


2.0 X 10"* 


Mg 


2.55 X 10"® 


P 


2.3 X 10"* 


CI 


1.4 X 10"* 



Table 3. Selected calculated fractional abundances 
at 10^ yr 



Species 






iVlOQei Z 




TT 

Xl 


o.D X lU 


-4 


0.0 X lU 


b 


H2 


0.50 




0.50 




CO 


8.9 X 10" 


-6 


9.2 X 10" 


6 





1.1 X 10" 


-5 


1.1 X 10" 


5 


C2 


3.7 X 10" 


-10 


4.4 X 10" 


10 


HCN 


6.1 X 10" 


-10 


7.5 X 10" 


10 


NH3 


4.2 X 10" 


-9 


4.3 X 10" 


9 


H2S 


1.8 X 10" 


-12 


2.0 X 10" 


12 


OH 


1.6 X 10" 


-9 


1.8 X 10" 


9 


C3H 


4.1 X 10" 


-11 


4.5 X 10" 


11 



designated Model 2, we removed all the surface reactions 
and ran the network again with the following additional 
changes: atomic hydrogen was not allowed to evaporate, 
and the accretion of H2 and He onto grains was halted, 
because we assume that these two species simply land on 
grains and quickly evaporate again. 

The initial fractional abundances for both models are 
given in Table [2l These are based on ele mental abundances 
considered bv IGarrod &: HerbstI (|2006[ ) based on unpub- 
lished work of Wakelam and Herbst. As compared with 
more standard low-metal elemental abundances, there is in 
general less depletion of the heavier elements. The constant 
density used is riH=2 x 10** cm"'^ and the cosmic ray ion- 



ization rate C is set to 1.3 x 10' 



-17 , 



We will use the same 



initial abundances and total density throughout the paper. 
Table [3] shows the comparison of fractional abundances of 
selected species for these two models at a time of 10^ yr. We 
can see that the difference is very limited; only the atomic 
hydrogen abundance is affected significantly; its abundance 
is greater in Model 1 because evaporation is allowed. But 
the extra atomic hydrogen does not make a significant dif- 
ference for other species. This agreement shows that the 
gas-phase abundances remain mostl y the same regardless 
of wh ether or not surface chemistry (jStantcheva fc Herbstj 
l2004f) is considered as long as the temperature is sufficiently 
low that evaporation is not efficient for heavy species and 
non-thermal desorption is not included. This latter con- 
straint is an important one because non-thermal desorption 
of surface ices may indeed be non-negligible for a few gas- 
phase species such as methanol. Nevertheless, in this first 
attempt to wed the Monte-Carlo microscopic approach with 
gas-phase chemistry, we ignore the process and focus atten- 
tion on the surface abundances. Under these conditions, we 
can infer that the details of surface kinetics are relatively 
insignificant for gas-phase abundances, whether or not they 
are handled by the rate-equation approach or by our more 
quantitative stochastic technique. The only important in- 
fluence of grains on gas-phase chemistry occurs by the ac- 
cretion of species onto grains, which eventually reduces the 
abundance of gas-phase species. 

4. Simulation IVIethod and Surface Models 

4.1. Calculation of Gas-Phase Chemistry 

Since the two different approaches to surface chemistry give 
very similar gas-phase abundances except for atomic hydro- 
gen, we can start with a rate-equation approach to the gas- 



phase abundances coupled with accretion of all species onto 
grains with unit sticking efficiency and no subsequent evap- 
oration of atomic hydrogen, as in Model 2. The assumption 
that ev aporation of H does not o ccur at all derives from th e 
work of IChang et all (|2005D and lCuppen fc HerbstI (|2005[) . 
who showed that at a grain temperature of 10 K the re- 
combination efficiency for molecular hydrogen is quite high 
for realistic surfaces. Nevertheless, it is not obviously true 
for all surfaces (see Table and iteration with the ac- 
tual Monte Carlo algorithm will be used to achieve better 
results. 

The error in the assumption that hydrogen atoms do not 
return to the gas phase after accretion can be calculated as 
follows. Suppose that Fevap{t) is the true time-dependent 
rate of atomic H evaporated from a grain. The evolution 
of the gas-phase abundance of atomic hydrogen Ug (H) can 
then be written as 



dngjH) 
dt 



^ formation^"* - ^ destruction^"'' 



^acc 

ng{H) + F, 



evap 1 



(1) 



where kacc is the accretion rate coefficient and the first two 
terms on the right-hand-side of the equation refer to gas- 
phase formation and destruction processes. Formal integra- 
tion of this equation yields 



-,{H)it) = J ^formation»""(/)dt' 

- JYI destruction^""'(t')dt' 



{l-p) k 

acc 

ng{H){t )dt 



where 



P 



(3) 



represents the ratio between the integrated evaporation rate 
and accretion rate of atomic hydrogen. This ratio can be 
easily determined during a Monte Carlo simulation, as dis- 
cussed later, by counting the number of evaporating and 
accreting atoms over the period of time. If p is much less 
than one, the error is small and the approximation of as- 
suming Fevap{t) to be zero is justified. If, however, p is 
large, the abundance of gas-phase atomic hydrogen can 
have a large error. Although this error will not affect other 
gas-phase abundances substantially, the incorrect hydrogen 
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abundance will not be suitable for the subsequent surface 
chemistry simulation, which depends on the accretion flux 
of H atoms. In this case, an iteration has to be performed to 
obtain the correct gas-phase H abundance as a function of 
time. As a first-order approximation, we assume the ratio 
of evaporating H to incoming H to be constant throughout 
the simulation and equal to p determined over the entire 
time scale of the simulation. This approach corresponds to 
taking the average evaporating effect into account. It is un- 
dertaken because the Gear routine used for the gas-grain 
calculation cannot solve ordinary differential equations with 
variables explicit in time. In other words, a variable l — p{t) 
factor cannot be inserted into the current gas-grain code. 
After p is obtained from the Monte Carlo simulation, we 
reduce the deposition rate coefficient of atomic hydrogen, 
kacc, to fcacc(l — p) and run the simulation again. If the 
new value of p is close to the previous one, the system is 
converged, otherwise the cycle is repeated with the new p 
until convergence is achieved. 

If p has no time dependence, the system converges to 
the true solution. However, p is typically time dependent 
and this time dependence cannot be included in the gas- 
grain network calculation. In order to gauge the resulting 
error, we define a time-dependent P2{t), which is the ra- 
tio of the number of evaporating H atoms to incoming H 
atoms counted during a "short" time period. We choose this 
"short" time to be 2 x 10^ yr, which is 1/100 of the total 
simulation time of 2 x 10^ yr. The time is short enough to 
determine the time dependence of p2 (t) and long enough to 
suppress most of the statistical noise. A shorter time scale 
(1 X 10'^ yr) shows some more noise while a longer time scale 
(5 X 10'^ yr) shows less. With the rate equation approach, 
this method is equivalent to the integration 



P2{t) 



t-2000 yr ^ ^^"-P 



(t )dt 



/t-2000 yr^accng{H){t )dt 



(4) 



If P2(t) is weakly dependent on time, t, then p has little time 
dependence, thus the resulting error is small. Otherwise, the 
error is large even if convergence has been achieved. With 
the current Gear algorithm used in the gas-grain model, we 
cannot use the parameter p2 {t) in the actual simulation. A 
new gas-grain algorithm would have to be designed to follow 
the gas-phase abundance of H with the time-dependent p2- 

4.2. Deposition of Gas-Phase Species 

Assuming the sticking coefficients for each species to be 1, 
the deposition (accretion) flux of species i onto each ab- 
sorption (lattice) site is equal to the product fiN, where fi 
is the incoming flux of species i in units of ML s~^, while 
N, the number of sites on the grain, is given by 



(5) 



where r is the radius of the grain and s is the surface site 
density. In our simulations, we start with bare olivine as 
the granular surface, and use a site density s of 1.5 x 10^^ 
cm~^, which is also used in our standard gas-grain network. 
The flux fi can be calculated by 



ng{i) SksTg 



TTrrii 



(6) 



where Tg is the gas phase temperature and rrii is the mass 
of species i. A factor of four, applicable to the case of spher- 
ical grains, is often used in the denominator. However, this 
factor is not used in the accretion term in our standard gas- 
grain network, and for consistency we do not use it here. 

If fi is independent of time, the deposition is an homo- 
geneous Poisson process, and the deposition interval At can 
then be calculated using 



At 



f^N ' 



(7) 



where X is a random number between and 1. In our 
simulations the deposition rate is time-dependent due to 
the time dependence of the gas-phase abundances of the 
impinging species. These gas-phase abundances are mean- 
field values, but the fiux is still stochastic in arrival time. 
In reality, there is already a distribution in the flux due to 
a distribution in velocity. The next deposition time, tnext, 
satisfies the equation. 



X = exp I - ^ "'"/,(t')iVdt'j , 



(8) 



with tprev the previous deposition time, and the process is 
referred to as an inhomogeneous Poisson one. It can easily 
be shown that Eq. ([8]) is also valid when fi is independent of 
time. With an analytical form of fi{t) known, tnext can be 
determined numerically. However, only a numerical solution 
is obtained from the gas-grain network calculation, thus 
interpolation is needed to approximate fi. The current gas- 
grain network outputs gas phase abundance of species i, at 
different time tj . We employ linear interpolation to get the 
abundance of species i at time t, which is between tj and 
tj+i. The flux fi{t) can then be calculated as 



(9) 



where A<j = tj+i ~ tj- We can formally write fi(t) for any 
t as, 

-S{t-t,+i)), (10) 



where S(t) is the step function and bi 
Combining Eqs. ^ and (fTO|) gives 



Atj 



X = exp(- 



{S{t ^tj)^S{t -tj+i))Ndt ). 
Let us call 

-t- {iprev J tnext) — j ^ ' ^ ^ ifiitj) bij{t ^j/)) 

"'tprei, j^Q 

{Sit ~t,)-Sit -t,+,))dt). 



(11) 



Then 
r = I{t 



prev 1 inext ) 7 



(12) 



(13) 



where r 



- \nX 
N ■ 



6 



Chang et al.: Stochastic gas-grain chemistry 



Eq. (1121) can be solved numerically in the following 
manner. Suppose tk < tprev < tk+i- First, we calculate 
I(tprevjtk+i)- If this is larger than r, we know that tnext 
must be less than tk+i- We can determine tnext by ex- 
pressing fi{t) in terms of tprev rather than in terms of tk- 
Integration then yields 



Table 4. Energies (K) used for different surfaces 



\l f i {tprev^ '^^^ik fi^prev) 



(14) 



where filtprev) = h{tk) + bikitprev - tk)- 

If I{tprev,tk+i) < r, then deposition must happen later 
than tk+i- We define a new quantity, r , 



r = r - I{tprev,tk+l) 



r — I{tk+l,tnext)- 



(15) 



(16) 



If integration from tk+i to tk+2 leads to a larger value 
than r , tnext < tk+2 and an expression for tnext similar 
to eq. (fT4|) can be obtained. If integration to tk+2 leads to 
a smaller value than r , tnext is larger than tk+2, and the 
sequence must be continued until tnext is bounded by two 
successive times tm and tm+i- 

To check for convergence, the original results of the gas- 
grain network are rerun twice with smaller time intervals, 
first 1/4 of the original Atj and then 1/4 of the second one. 
Convergence is achieved in all cases. 

4.3. Surface Models 

A detailed de scription of fiat and rough surface models has 
been given in ICuppen fc HerbstI ^05). Here we start with 
either a fiat surface or their surface D, which is the roughest 
one and is generated by another Monte Carlo program. The 
diffusion barrier. Eh, and the desorption energy, E]j, for the 
two diffusing species - atomic hydrogen and atomic oxygen 
- are in our models dependent on the topological rough- 
ness. Starting from either a fiat or rough bare olivine sur- 
face, a mantle gradually develops. The formation of surface 
species, especially water ice, on the surface has two conse- 
quences. First, the chemical nature of the surface changes, 
resulting in different energetics and therefore different val- 
ues for Eb and Ed- Secondly, the surface topology changes, 
as more and more species are formed and roughness devel- 
ops even from the initially flat surface due to fluctuations. 
Thus, Eii and Eo are inherently stochastic in the simula- 
tion in that they change for the site at the top of a given 
column as the stochastic calculation proceeds. For simplic- 
ity, we assume that only the olivine substrate and the H2O 
molecules (the main mantle constituent) influence E^ and 
E]j, either by lateral interactions or by a vertical interac- 
tion with the substrate. 

The parameters Eh and Ejj are given by 



Eh = El + iEl+]E_ 



Ed — Ed 



lEl 



3 Eh 



(17) 
(18) 



respectively. Here E^ and E^ are the diffusion barrier 
and desorption energy on a surface without lateral bonds, 
while E]^ and E\ are the lateral bonds imposed by hor- 
izontal olivine and H2O neighbors respectively. In this 
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simulation. El has a value of Q-IEd (jCuppen fc Herbs^ 
2005). The indices i and j are the number of olivine and 
H2O neighbors respectively. Since species other than H 
and O have large diffusion barriers even without lateral 
bonds (|Ruffie fc Herbstl[2000l) . we ignore their movement in 
the temperature range studied here (10-15 K). The evapo- 
ration of O is ignored. 

Table m summarizes the parameters used for H and O on 
the ice and olivine surfaces. The diffusi on barrier and des - 
orption energy of H on olivine are from lKatz eFaLl (|1999[ ). 
The desorption energ y for H on water ice is not well con- 
strained. HoUenbach fc Salpeterl (jl970[ ) found a value of 450 
K, Buch fc Czerminskil (|1991[ ) found a value of 500 K, while 
Al-Halabi et all ()2002[) found a value of 400 K. We take 
the av erage of these measurements. Recently, iPerets et al.l 
()2005D deduced a much higher value of 7 20 K. The oxygen 
values are from iTielens fc HagenI ()1982f ). For comparison, 
we also simulate a second system by ignoring the lateral 
bonds and only considering the change of the chemical na- 
ture of the surface, i.e. the vertical bond. This system starts 
from a flat surface. Although it gains some roughness from 
fluctuations as the ice layers grow, the system is kinetically 
equivalent to a flat system since there is no penalty for 
going uphill and no gain in falling into a valley. 

4.4. Surface reactions with activation energy barriers 

In gas-grain codes based on a rate equation or master equa- 
tion approach to the surface chemistry, the common method 
to calculate diffusive rate coefficients for reactions with an 
activation energy barrier is to multiply the diffusion rate 
by a tunneling factor: 



2a 



= exp ( - — y^lEJiJ] 



(19) 



where a is the width of the rectangular barrier, typically 
chosen to be 1 A, /i is the reduced mass, and E g is the re- 
action activation energy (|Hasegawa et al.lll992f) . Note that 
the word "diffusion" is normally used to refer to motion 
over a whole grain, whereas the word "hopping" refers to 
motion from one potential minimum to the next. A physical 
interpretation of the rate coefficient is that the two species 
have a probability equal to the tunneling expression to re- 
act with each other if they diffuse into the same potential 
minimum. Thus, the rate coefficients are linearly dependent 
on the diffusion rate. 

Unlike the typical situation in the gas phase, two species 
that meet each other on a surface can stay in each other's 
vicinity for quite a while depending on the hopping rate to 
nearby lattice sites. During this time they can have many 
chances to react. It is therefore more physical to model re- 
actions with barriers as processes competitive with hopping 
out of the potential minimum. Such an approach has been 
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done by lAwad et al.l (|2005f ) using rate equations. Here we 
incorporate this approach by considering the competition 
among Poisson processes. 

In particular, the competition among reaction, hopping 
or evaporation can be implemented as follows. The reaction 
rate with tunneling can be written as 



i/i exp 



2a 



(20) 

where vi is the attempt frequency for reaction. The com- 
peting rates for hopping and evaporation of the two species 
are given by 

bi = V2exp (-Eb/T) , 

and 

62 = i^scxp {-Ed/T) , 



(22) 



respectively, where T is the surface temperature, 1^2 is the 
attempt frequency for hopping, and is the attempt fre- 
quency for evaporation, and the diffusion barrier and evap- 
oration energy pertain to either of the two species. The 
stochastic algorithm for diffusion treats all directions with 
equal probability. Note that hopping is treated classically 
rather than by quantum mechanical tunneling because, un- 
like reaction, hopping is thought to occur over a broad and 
shallow potential. Quantum chemical calculations are be- 
ing undertaken to verify this point (Woon, private com- 
munication). Assuming that only one species is moving, 
the probability for reaction to occur instead of diffusion or 
evaporation is 



P 



(23) 



We assume 1^1 = 1^2 = I's = 10^^ s^^, but since there three 
different types of processes the attempt frequencies can be 
different. The competition is implemented in the Monte 
Carlo algorithm as follows: a random number between 
and 1 is generated. If this number is smaller than p the reac- 
tion occurs. Otherwise evaporating or hopping takes place. 
A competition between these two processes must also be 
evaluated. 



4.5. Monte Carlo Simulation 

The Monte Carlo algorithm starts with determining the 
first deposition time for H, O, and CO, as explained in 
Section [42l Then the first deposition is executed by plac- 
ing the particle on a randomly picked site of the lattice, 
which is chosen to have a lattice of size 100 x 100. This 
size corresponds to a rather small grain of 0.015 /zm, and 
a mean distance between atoms in adjacent lattice sites of 
?a 5 A. Typical activation energy barriers are shorter than 
this in width; the widths of the diffusion barriers are cur- 
rently under investigation (Woon, private communication). 
We have already shown that the lattice size used is large 
enough that there is little size depend ence to the results 
for both flat and rough olivine surfaces (jChang et al.ll2005t 
ICuppen fc Herbstll2005[ ). Nevertheless, we redid our simu- 
lations for lattice sizes of 50 x 50 and 200 x 200 lattices at 
15 K to confirm this point. For much smaller grains, the so- 
called accretion limit is reached, where the average mantle 
abundances of reactive species is less than unity, leading to 



a drop in reaction efficiency. For everything but the surface 
stochastic simulations, the dust-to-gas number density is 
taken to be 1.32 x lO^^^ and the grain size the nominal one 
of radius O.l/im. 

As the clock moves forward, other species are deposited 
on the lattice, and species already on the lattice can hop 
from site to site or evaporate. By hopping on a rough sur- 
face, we mean motion on the top level of any occupied site, 
whether is is occupied by the original silicate or by any 
other species. Thus hopping can include climbing or falling 
in a vertical sense. For very rough surfaces, the climbing 
or falling length can be up to 20-30 molecules. Evaporation 
and hopping are assumed to be Poisson distributed, with 
time intervals between events At = — ^'^^^"^ , where A" is a 
random number between and 1, and b is the occurring 
rate of that process as given by Eqs. pTjl and (|22p. For 
H atoms, which are allowed to both hop and evaporate, 
these two processes are treated first as one combined pro- 
cess with a rate 61 -I- &2 and then another random number 
X is called to decide which one will occur first. If X is 
smaller than ^^^^^ , the hydrogen atom hops and it evapo- 
rates in the other case. The total number of H atoms that 
evaporate is counted to calculate p. Oxygen atoms are only 
allowed to hop. The time interval between hopping is solely 
determined by the hopping rate. 

If a hydrogen or oxygen atom hops to a site where there 
is already one species that can react with it without a bar- 
rier, the reaction happens immediately. If either atom hops 
into the same site as a species with which it can react with 
an activation energy barrier, the competition among reac- 
tion and other processes occurs as discussed in the previ- 
ous section. Reaction products other than II2 remain on 
the surface and can cover species on the same sites. Some 
reaction products can react further with species they cover. 
For instance, if OH is produced by a reaction in which H 
hops into a site occupied by an O atom, and there is an 
H2CO molecule lying under the O atom, then another re- 
action is allowed to take place, producing HCO and H2O 
(see Table [J). 

If a gas-phase species lands atop a bound surface 
species, reaction can occur according to the Eley-Rideal 
mechanism. No distinction is made between reaction by 
hopping (Langmuir-Hinshelwood mechanism) and by land- 
ing. If reaction does not occur upon landing, the incoming 
particle covers the one already present. No reaction is al- 
lowed between species covered by at least one other species. 
For instance, assume that an oxygen atom is buried by a 
CO molecule. After being buried, the oxygen atom cannot 
react with H or O even if they occupy the site above the 
CO. 

The Monte Carlo simulation is performed by executing 
the following sequence of steps: 

1) The gas-grain code is run for Model 2 through 10^ yr, 

2) The gas- phase results for H, O, CO are used to run the 
Monte Carlo simulation of deposition, surface chemistry, 
and evaporation of atomic hydrogen for 2 x 10^ yr, 

3) An average value of p over the full time of the stochastic 
calculation is obtained, 

4) If p is close to zero or is converged, the calculation is 
stopped, otherwise: 

5) The gas-grain code is rerun with a deposition rate 
reduced by a factor of 1 - p, and the simulation is also run 
with a reduced deposition rate for H, 
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Fig. 2. p2{t) vs time t for a surface-adsorbate system with- 
out lateral bonds at 15 K. The solid line is for p2 after 
Fig. 1. Fractional abundances of gas-phase species that ac- iteration, while the dotted hue is for p2 calculated in the 
Crete onto grains for Model 2. Solid lines are for 10 K while first Monte Carlo simulation, 
dotted lines are for 15 K. 



6) The procedure then returns to step 3. 

Note that the gas-grain code is used to obtain gas- 
phase results; surface results are subsequently obtained by 
stochastic simulation. The Monte Carlo code is run for a 
shorter time than the gas-grain code for two main reasons: 
(i) it is very time-consuming computationally, and (ii) the 
parameter p2 {t) becomes significantly time-dependent after 
2 X 10^ yr. Finally, the iterative procedure used is forced 
on us by the current structure of the gas-grain code, which 
does not allow parameters with explicit time dependence. 



5. Results and Discussion 

The simulation starts in earnest with runs at 10 K and 15 
K in which the fluxes of hydrogen atoms, oxygen atoms, 
and CO striking the grain as a function of time are taken 
from Model 2. Figure [1] shows the fractional abundances of 
gaseous H, O and CO as functions of time at 10 K and 15 K 
obtained from the Model 2 run. As can be seen, the atomic 
hydrogen abundances do not change much over the first 
10^ yr, while the decline in the atomic oxygen abundance 
turns on after 10^ yr. A significant amount of gas-phase 
CO is produced after 10'^ yr. These abundances are those 
used to calculate fluxes for the Monte Carlo simulation of 
the surface chemistry. In this run, both surface chemistry 
and evaporation of atomic hydrogen occur. As discussed in 
Section [4Tl the evaporation of atomic hydrogen is followed 
stochastically and averaged over time to determine the pa- 
rameter p. If this parameter is significantly different from 
zero, then further runs are necessary with a reduced hydro- 
gen flux to achieve convergence. The closer the value of p 
to unity, the slower the convergence. 

We first simulate the surface chemistry starting with a 
flat olivine surface and using a surface-adsorbate system 
with no lateral bonds as the ice develops. We then redo 
the simulation starting from a rough olivine surface and 
containing lateral bonds. The latter surface, which is rough 
at all stages, will just be referred to as "rough" . 



5.1. Surface without lateral bonds 

At 10 K, p is initially calculated to be 3 x 10^^, which 
is significantly less than 1. Thus, the gas-phase fractional 
abundances shown in Fig. [T] are also reasonably good solu- 
tions of the gas-grain network with evaporation. However, 
at 15 K, p is 0.37, which is comparable to 1, because evap- 
oration is more important. Following the prescription in 
Section [4Tl we reduce the hydrogen-atom accretion rate to 
0.63 of the previous value, and perform the gas-grain simu- 
lation a second time with the Monte Carlo approach to the 
surface chemistry. In this simulation, p now is calculated 
to be 0.40, which indicates that reasonable average con- 
vergence has been reached. Moreover, we can look at the 
time-dependent p2 (t) to check that it is reasonably indepen- 
dent of time, another test of convergence. Figure [H shows 
P2(t) as a function of time for 15 K. From this figure, we see 
that the error of our approach is small since p2 (t) converges 
and does not change much with time, although there is an 
increase after 10^ yr. The increase of P2{t) derives from the 
decrease of O in the gas phase, which decreases its flux onto 
grain surfaces, so that more H atoms can evaporate before 
combining with surface O. 

The gaseous fractional abundance of atomic hydrogen 
after the initial Monte Carlo simulation and its iteration 
at 15 K is shown in Fig. [31 As expected, the H frac- 
tional abundance increases after the initial surface simu- 
lation since evaporation is allowed. For example, at 10*^ 
yr, the H fractional abundance exceeds 10"'' whereas it is 
closer to 7 X 10~^ for Model 2. After the iteration, however, 
the atomic H abundance has become closer to its original 
abundance at all times. Even with the increased gas-phase 
H abundance after the initial simulation, the abundance for 
each surface species changes less than 15% at 15 K. This in- 
dicates that the change in the gas-phase atomic hydrogen 
abundance due to the change in H evaporation from the 
grain surface has little effect on the surface abundances. 
The other gas-phase abundances do not change much ei- 
ther because atomic hydrogen only plays a minor role in 
the gas-phase chemistry. 

Figure m shows the surface abundances of CO, H2CO, 
CH3OH and CO2 on grain surfaces as functions of time at 
10 K and 15 K. After some initial fluctuations, especially 
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Table 5. Mantle populations (ML) 2 x 10^ yr (no lateral bonds). 



1000 
time(yrs) 



Fig. 3. Fractional abundance of gas-phase atomic hydrogen 
at 15 K. The solid line refers to the first simulation with 
the Monte Carlo method while the dotted line refers to a 
subsequent iteration. 
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Fig. 4. Mantle abundances for selected species in monolay- 
ers as a function of time on a surface-adsorbate system with 
no lateral bonds. Circles are for CO, triangles pointed up- 
ward are for CO2, squares are for H2CO, while diamonds 
are for CH3OH. The left panel refers to results at 10 K 
while the right panel refers to results at 15 K. 



at 10 K, the population of each carbon-bearing species (in 
monolayers per grain) increases monotonically with time. 
At 10 K, methanol exceeds CO after 500 yr and eventually 
becomes the most abundant carbon-bearing surface species, 
while CO2 is also abundantly produced. At 15 K, on the 
other hand, methanol is always the least abundant of the 
four carbon-bearing species shown. 

Table [5] shows the populations of all the mantle species 
at a time of 2 x 10^ yr. It is easily seen that the mantles are 
dominated by water ice, with molecular oxygen second in 
abundance. The total number of species at 15 K is slightly 
larger because of the larger accretion rate at higher temper- 
atures and the still minimal evaporation rate. As already 
seen in Fig. 21 the abundance of methanol is highly depen- 
dent on temperature; its abundance drops by a factor of 
seven as the temperature increases from 10 K to 15 K. The 
higher temperature makes the formation of methanol less 
efficient because diffusion rates increase quickly at higher 
temperatures while the tunneling reaction rates for the crit- 
ical reactions with activation energy - H -I- CO and H -|- 
H2CO - are independent of temperature. Thus, competi- 
tion favors hopping of H over reaction as the temperature 
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Fig. 5. Abundance per monolayer in fractions of a mono- 
layer for selected species at each layer on a surface without 
lateral bonds at 2 x 10^ yr. Monolayers are numbered from 
zero at the initial surface. The symbols representing the 
species are the same as in Fig. S) The left and right panels 
represent results at 10 K and 15 K, respectively. 



increases. This decrease in the methanol abundance with 
increasing temperature is not exactly what is measured 
in the laboratory, although the conditi ons there are some - 
what different from those studied here (|Fuchs et al.ll2007t ). 
Our analysis for the calculated temperature dependence of 
methanol also leads to more HCO and H3CO at 10 K than 
at 15 K. In addition, the increasing temperature makes hy- 
drogen atoms hop faster, thus the activation energy-less 
reactions with HCO and H3CO become more efficient. For 
HCO and H3CO, there is another mechanism for destruc- 
tion: O atoms become mobile at 15 K, and can react with 
both radicals. Unlike other simulations, much CO2 is pro- 
duced in our calculation, mainly because products from one 
reaction can react further with species under them; thus, 
species that are not mobile can react with each other. In 
this case, the reaction of importance is that between CO 
and OH, the barrier of which is low enough that the prob- 
ability of reaction is near unity. For the case of O2, its for- 
mation by surface recombination of oxygen atoms is helped 
by an increase in temperature since the surface atoms move 
around significantly more quickly at 15 K. At 10 K, there is 
a significant amount of atomic oxygen in the grain mantles, 
most of which is buried and not able to react further. 

Figure [5] shows the abundances of the major carbon- 
bearing species in individual monolayers, starting from the 
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Fig. 6. Mantle abundances for selected species in monolay- 
ers as a function of time on a rough surface. Circles are for 
CO, triangles pointed upward are for CO2, squares are for 
H2CO, while diamonds are for CH3OH. The left panel refers 
to results at 10 K while the right panel refers to results at 
15 K. 



surface and moving upward. These differential abundances 
of each carbon-containing species typically increase initially 
as a function of height from the original surface but peak 
before the top layers are reached at both 10 K and 15 K. 
The bottom layers hardly have any molecules with carbon, 
because gas-phase CO takes some time to produce. The top 
layers contain fewer species overall because these layers are 
not yet fully occupied. 

5.2. Rough surface 

The calculated p at 10 K is about 1 x 10^^, which is much 
smaller than unity, so that no iteration is required. At 15 
K, p is 0.14, which is still smaller than 1. But we perform 
an iteration anyway. After one iteration, p is still 0.14, so 
the calculation is converged. 

The time-dependent surface abundances of CO, H2CO, 
CH3OH and CO2 are shown in Fig. El At 10 K, the figure 
is similar to Fig. 31 which represents the analogous results 
on a surface-adsorbate system without lateral bonds, while 
at 15 K, on the other hand, more methanol is produced 
than in the system without lateral bonds, because there are 
stronger binding sites due to lateral bond interactions and 
these stronger binding sites help to add hydrogen atoms to 
molecules when there is a reaction barrier, since the com- 
petitive diffusion rate is lowered. 

The surface abundances of all studied species at 2 x 10^ 
yr are shown in Table [51 Compared with Table [5l there 
is not much difference at 10 K except that more hydro- 
gen atoms are trapped on the surface. Indeed, for any 
reasonably-sized grain, the mantle population of 0.019 cor- 
responds to far greater than 1 atom per grain. At 15 K, 
however, as in Fig. [6l one can see that more methanol is 
produced than without lateral bonds. 

The abundance per monolayer of the carbon-bearing 
species at different layers is shown in Fig. [7] for a time of 
2 X 10^ yr. For 10 K, this gives a very similar result to 
the system without lateral bonds in Fig. [51 At 15 K, more 
methanol is formed at each layer because of the slower hop- 
ping rates caused by the lateral bond interaction. 



Table 6. Mantle populations (ML) 2 x 10^ yr (rough surface) 
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Fig. 7. Abundance per monolayer in fractions of a mono- 
layer for selected species at each layer on a rough grain 
surface at 2 x 10^ yr. Monolayers are numbered from zero 
at the initial surface. The symbols representing the species 
are the same as in Fig. 5. The left and right panels repre- 
sent results at 10 K and 15 K, respectively. 



At the rough surface at 15 K, the populations of species 
as a function of monolayer have a distinct feature, in that 
they drop to zero abruptly, unlike the other cases. This 
effect happens because the oxygen atoms are able to hop 
efficiently out of sites with small lateral bond interaction, 
located above surrounding sites, to sites with larger lat- 
eral bond interaction, located in valleys below surrounding 
sites, where they are trapped. So, sites that are lower in ver- 
tical height above the original surface tend to have a larger 
probability to be occupied by oxygen atoms. These atoms 
can remain or react with other atoms to fill the up the val- 
leys. Thus, the diffusion-reaction mechanism automatically 
"smooths" the surface strongly at 15 K. 

In addition to the rapid dropping in monolayer abun- 
dance of the four major carbon-bearing species, there is a 
second distinct feature at 15 K. Both the monolayer abun- 
dances of CO and II2CO increase strongly towards the high- 
est monolayers before suddenly dropping. This effect occurs 
because it is relatively difficult to bury CO and II2CO un- 
der higher layers due to their reactivity with OH radicals. 

5.3. Comparison with observations 

Our primary interest in this work has been to introduce 
a new method to combine rate equations and microscopic 
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Monte Carlo simulations in a gas-grain model of the chem- 
istry of cold interstellar cores. Toward this goal, a small 
grain-surface chemical network was chosen for simplicity. 
However, we would still like to compare our results with as- 
tronomical observations. Since, again for simplicity, we do 
not consider non-themal desorption mechanisms here, as we 
have done recent ly with our latest gas-grain rate-equation 
code ()Garrod et a l. 2006 ). it is unnecessary to compare our 
gas-phase results with o bservations since the r esults are 
very similar to those of iRufHe fc HerbstI (|2000f ). Rather, 
it makes more sense at this stage to compare our mantle 
results with infrared absorption measurements. These re- 
sults are insensitive to non-thermal desorption, which only 
affects mantle abundances slightly. We are particularly in- 
terested in determining whether or not a rough surface 
model can help to obtain abundances closer to the observa- 
tions. In Table [7] we list calculated and observed fractional 
abundances for four ices: H2O, CO2, CO, and CH3OH. The 
water abundance is tabulated as the fractional abundance 
with respect to nn while the abundances of the other ices 
are listed as percentages of the water ice abundance. Our 
theoretical results are calculated at 2 x 10^ yr. For compar- 
ison, master equat ion results at 10 K for a tim e of 3.2 x 10^ 
yr are also listed (jStantcheva fc HerbstI 2004 ). The o bser- 
vational da t a for W 33 A are from lDartois et al.l (|1999D and 
l2000f) . while the data for Elias 16 are from 
2OOOD . 

As compared with the master equation results, which 
are based on a macroscopic stochastic approach in which 
no distinction is made for local site differences and differ- 
ences according to monolayer, the Monte Carlo results are 
a distinct improvement. The master equation abundance 
for CO is very small at 10 K, although better results are 
obtained at 20 K. As for CO2, the results are low at both 
10 K and 20 K. For the Monte Carlo models, agreement is 
reasonably good regardless of the type of surface model ex- 
cept perhaps for CO. The agreement for CO2 is especially 
important given the ubiquity and high abundance of this 
ice component. The W33A results are closer to our 10 K 
model predictions whereas the Elias 16 results are closer 
to the 15 K predictions. The higher temperature leads to 
more CO and less methanol. Still, CO is somewhat under- 
produced in the models. For the rough surface, increasing 
the lateral bond strength to OAEd from its standard value 
of O.lii^D can increase CO from 0.75 percent to 3.0 percent 
of the H2O ice while keeping other species in good agree- 
ment with observation at 10 K. 



6. Conclusions 

For the first time, a gas-grain chemical simulation of a dense 
interstellar cloud core has been performed with a combined 
rate equation-microscopic Monte Carlo approach. In this 
approach, the local characteristics of the surface and all 
monolayers accreted onto it can be followed. The surface 
network is chosen to be the H, O and CO system, which 
has been widely studied and is the key network to form 
methanol on grain surfaces. Our results for the total mantle 
abundances of the carbon-bearing species CO, CO2, and 
CH3OH as well as for water are in reasonable agreement 
with observation in two well-studied protostellar sources. 
In fact, the agreement is superior to that achieved by use of 
a macroscopic Mont e Carlo approach known as the master 
equation treatment (jStantcheva fc Herbstll20d^ . 



The approach utilized here is based on a number of sim- 
plifications made possible mainly by the low temperatures 
(10 K, 15 K), the starting chemical abundances, and an un- 
usual aspect of the gas-phase chemistry. First, if we start 
with hydrogen mainly in molecular form, we need not con- 
sider subsequent H2 formation since it is a small effect. 
Secondly, at the low temperatures used here, the only reac- 
tive species that can evaporate from grain surfaces is atomic 
hydrogen. Thirdly, the residual atomic hydrogen abundance 
in the gas does not strongly affect the abundances of all 
other gas-phase species. These simplifications allow an it- 
erative solution in which gas-phase abundances computed 
initially do not change significantly except for the case of 
atomic hydrogen. Our iterative procedure involves starting 
with a gas-grain calculation in which no surface chemistry 
or H evaporation is allowed. The time-dependent abun- 
dances of H, O, and CO arc then allowed to deplete onto 
grain surfaces and the surface chemistry and evaporation of 
H are treated by the continuous-time random walk Monte 
Carlo procedure. The evaporation rate of hydrogen is used 
to refine the abundance of atomic hydrogen in the gas by 
an average correction over the time period of the calcu- 
lation, an approximation made possible by the neglect of 
H2 formation. This refined abundance is then used, if nec- 
essary, to redo the Monte Carlo calculation of the surface 
chemistry until convergence is achieved. If the evaporation 
of atomic hydrogen from grain surfaces is strongly time de- 
pendent, a more complex bootstrap procedure would be 
called for. Were other species to evaporate, as occurs at 
higher temperatures, and were non-thermal desorption to 
be included, a process that occurs even at low tempera- 
tures, we would probably have to devise a more intricate, 
non-iterative, method of integrating the gas-phase abun- 
dances simultaneously with the Monte Carlo determination 
of surface abundances. Such improvements arc under active 
consideration, including a n approach similar to the effective 
rate coefficient method (jCarrod et al.l 120061 : IChang et al.l 
l 2006f ). Finally, only a small number of surface processes 
are currently used; inclusion of a much larger number of 
processes by the Monte Carlo procedure would be exceed- 
ingly time-consuming whether or not changes in gas-phase 
abundances are also calculated. Whether a suitable proce- 
dure that will allow us to remove all of these simplifications 
can be achieved in the near future remains uncertain at his 
point. 
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